#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
mat identity_arma(int iNN) {
// Create & return an identity matrix
  mat AA = eye( iNN, iNN ) ;
  return AA ;
}


// [[Rcpp::export]]
mat invert_cpp(colvec expmeanval, colvec share, mat expmu, mat oo, sp_mat sharesum, uvec marketForProducts) {

int ii = 0;
double norm_max = 1;
double tol_inner = 1.e-14;
mat expmeanval0 = expmeanval;
mat numer = (expmeanval0*oo)%expmu;  
mat sum1 = sharesum*numer; 
mat sum11 = 1/(1+sum1);  
mat denom1 = sum11.rows(marketForProducts);
mat SS = numer%denom1; 
colvec PROB = mean(SS,1);  
colvec expmeanval1 = expmeanval0%share/PROB;
colvec t = abs(expmeanval1-expmeanval0);


while ((norm_max > tol_inner) && (ii < 2500)) {
  
numer = (expmeanval0*oo)%expmu;  
sum1 = sharesum*numer; 
sum11 = 1/(1+sum1);   
denom1 = sum11.rows(marketForProducts);
SS = numer%denom1; 
PROB = mean(SS,1);  


expmeanval1 = expmeanval0%share/PROB;
t = abs(expmeanval1-expmeanval0);
norm_max = as_scalar(max(t,0));
expmeanval0 = expmeanval1;
++ii;
}
return expmeanval0;
}


// [[Rcpp::export]]
List ind( mat expmeanval, mat expmu, mat oo, sp_mat sharesum) {
  
mat numer=(expmeanval*oo)%expmu;
mat sum1=sharesum*numer;
mat sum11=1/(1+sum1);
return List::create( Named("sum11") = sum11, Named("numer") = numer ) ;
}  


// [[Rcpp::export]]
mat objective( mat x, double K, double numProdsTotal, mat W) {
  
mat g = x.rows((2*K+2+numProdsTotal),x.n_rows-1);  
mat f = g.t()*W*g;

return f ;
}  



// [[Rcpp::export]]
colvec hess0_cpp(mat simShare, int index_start, int index_end, int rr) {
  
colvec simS = simShare.submat(index_start-1, rr-1, index_end-1, rr-1);

return simS ;
}

// [[Rcpp::export]]
mat hess01_cpp(mat simS_mat) {
  
int rows = simS_mat.n_rows;  
mat simS_diag = eye(rows,rows)%repmat(simS_mat,1,rows);

return simS_diag ;
}



// [[Rcpp::export]]
double hess_cpp( mat multip, mat simS) {
  
double sumprod_mpsimS = as_scalar(trans(multip)*simS);

return sumprod_mpsimS ;
}  


// [[Rcpp::export]]
mat hess2_cpp(mat x_char_index, mat v_rr, mat ooo1, mat simS_mat) {
        
mat xsimSx = x_char_index - (ooo1*(trans(simS_mat)*x_char_index));
mat xsimSxv = xsimSx%(ooo1*trans(v_rr));

return xsimSxv ;
} 

// [[Rcpp::export]]
mat hess3_cpp(mat xsimSxv, mat ooo, mat simS_mat) {
        
mat dSdtheta2rr = (simS_mat*ooo)%xsimSxv;

return dSdtheta2rr ;
} 



// [[Rcpp::export]]
mat hess4_cpp(mat simS_mat, mat multip, mat dSdtheta2rr, double sumprod_mpsimS, mat ooo, int nn, mat dL2ddeltadtheta_old) {
        
mat dL2ddeltadthetarr = -simS_mat*trans(multip)*dSdtheta2rr - sumprod_mpsimS*dSdtheta2rr + (multip*ooo)%dSdtheta2rr;   
mat dL2ddeltadtheta_new = dL2ddeltadtheta_old + dL2ddeltadthetarr/nn;

return dL2ddeltadtheta_new ;
} 


// [[Rcpp::export]]
mat hess5_cpp(mat dL2ddeltadtheta_old, mat dL2ddeltadthetarr, int nn) {
        
mat dL2ddeltadtheta_new = dL2ddeltadtheta_old + dL2ddeltadthetarr/nn;

return dL2ddeltadtheta_new ;
} 


// [[Rcpp::export]]
mat hess6_cpp(mat multip, mat ooo, mat dSdtheta2rr, mat xsimSxv, double sumprod_mpsimS, mat x_char_index, mat v_rr, mat matrix_ones, int nn, mat dL2dtheta22_old) {
        
mat dL2dtheta22rr = trans((multip*ooo)%dSdtheta2rr)*xsimSxv + sumprod_mpsimS*(-trans(dSdtheta2rr)*x_char_index%(matrix_ones*trans(v_rr)));
mat dL2dtheta22_new = dL2dtheta22_old + dL2dtheta22rr/nn;

return dL2dtheta22_new ;
} 


// [[Rcpp::export]]
mat hess_pre2_end_cpp(mat dL2ddeltadtheta2_old, mat simS_mat, mat multip, mat ooo, mat dSdtheta2rr_mat, mat xsimSxv_mat, double sumprod_mpsimS, mat x_char, int index_start, int index_end, mat v, int K, int nD, int rr, int nn) { 
  
mat dL2ddeltadtheta2rr ;   
mat dL2ddeltadtheta2_new = dL2ddeltadtheta2_old ;
  
for ( int dd = 1; dd < (nD+1); ++dd ) {  
dL2ddeltadtheta2rr = -simS_mat*trans(multip)*dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1) - sumprod_mpsimS*dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1) + (multip*ooo)%dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1) ;      
dL2ddeltadtheta2_new.submat(index_start-1,((dd-1)*(K+1)), index_end-1, (dd*(K+1)-1)) = dL2ddeltadtheta2_old.submat(index_start-1,((dd-1)*(K+1)), index_end-1, (dd*(K+1)-1)) + dL2ddeltadtheta2rr/nn ;    
} 
 
return dL2ddeltadtheta2_new ;
} 


// [[Rcpp::export]]
mat  hess_pre_end_cpp(mat dL2dtheta2rr, mat dL2dtheta2_old, mat multip, mat ooo, mat dSdtheta2rr_mat, mat xsimSxv_mat, double sumprod_mpsimS, mat x_char, int index_start, int index_end, mat v, int K, int nD, int rr, int nn) { 
  
mat dL2dtheta2_new = dL2dtheta2_old ;
  
for ( int dd = 1; dd < (nD+1); ++dd ) {  
dL2dtheta2rr.submat(((dd-1)*(K+1)), ((dd-1)*(K+1)), (dd*(K+1)-1), (dd*(K+1)-1)) = trans(((multip*ooo)%dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1)))*xsimSxv_mat.cols((K+1)*(dd-1), dd*(K+1)-1) + sumprod_mpsimS*(-trans(dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1))*x_char.rows(index_start-1, index_end-1)%(ones(K+1,1)*trans(v.submat((K+1)*(dd-1), rr-1, dd*(K+1)-1, rr-1)))) ;
dL2dtheta2_new.submat(((dd-1)*(K+1)), ((dd-1)*(K+1)), (dd*(K+1))-1, (dd*(K+1))-1)= dL2dtheta2_old.submat(((dd-1)*(K+1)), ((dd-1)*(K+1)), (dd*(K+1)-1), (dd*(K+1))-1)  + dL2dtheta2rr.submat((dd-1)*(K+1), ((dd-1)*(K+1)), (dd*(K+1)-1), (dd*(K+1)-1))/nn ;
}
    
return dL2dtheta2_new ;
} 

// [[Rcpp::export]]
mat  hess_end_cpp(mat dL2dtheta2_old, mat multip, mat ooo, mat dSdtheta2rr_mat, mat xsimSxv_mat, double sumprod_mpsimS, mat x_char, int index_start, int index_end, mat v, int K, int nD, int rr, int nn) { 
  
mat dL2dtheta2rr ; 
mat dL2dtheta2_new = dL2dtheta2_old ;
  
for ( int dd = 1; dd < (nD+1); ++dd ) {  
for ( int ddd = 1; ddd < (nD+1); ++ddd ) {  
if (ddd!=dd) {   
dL2dtheta2rr = trans((multip*ooo)%dSdtheta2rr_mat.cols((K+1)*(ddd-1), ddd*(K+1)-1))*xsimSxv_mat.cols((K+1)*(dd-1), dd*(K+1)-1) + sumprod_mpsimS*(-trans(dSdtheta2rr_mat.cols((K+1)*(ddd-1), ddd*(K+1)-1))*x_char.rows(index_start-1, index_end-1)%(ones(K+1,1)*trans(v.submat((K+1)*(dd-1), rr-1, dd*(K+1)-1, rr-1)))) ;   
dL2dtheta2_new.submat(((dd-1)*(K+1)), ((ddd-1)*(K+1)), (dd*(K+1))-1, (ddd*(K+1))-1) = dL2dtheta2_old.submat(((dd-1)*(K+1)), ((ddd-1)*(K+1)), (dd*(K+1))-1,(ddd*(K+1))-1) + dL2dtheta2rr/nn ;
}
}
}

return dL2dtheta2_new ;
} 
       

// [[Rcpp::export]]
mat hess_endend_cpp(mat hess_old, mat dL2dtheta2, mat dL2ddeltadtheta2, mat dL2ddelta2, int nD, int K, int numProdsTotal) { 
  
mat hess = hess_old ;  
 
// HESSIAN FOR THETA2, THETA3, THETA4 
for ( int dd = 1; dd < (nD+1); ++dd ) {
hess.submat(dd*(K+1), dd*(K+1), (dd+1)*(K+1)-1, (dd+1)*(K+1)-1) = dL2dtheta2.submat((dd-1)*(K+1), (dd-1)*(K+1), dd*(K+1)-1, dd*(K+1)-1) ;    
} 
 
for ( int dd = 1; dd < (nD+1); ++dd ) {  
for ( int ddd = 1; ddd < (nD+1); ++ddd ) {  
if (ddd!=dd) {     
hess.submat(dd*(K+1), ddd*(K+1), (dd+1)*(K+1)-1, (ddd+1)*(K+1)-1) = trans(dL2dtheta2.submat((dd-1)*(K+1), (ddd-1)*(K+1), dd*(K+1)-1, ddd*(K+1)-1)) ;    
}
}
} 
 
// HESSIAN FOR DELTA AND THETA2dd
for ( int dd = 1; dd < (nD+1); ++dd ) {
hess.submat((nD+1)*(K+1), dd*(K+1), (nD+1)*(K+1)+numProdsTotal-1, (dd+1)*(K+1)-1) = dL2ddeltadtheta2.cols((dd-1)*(K+1), dd*(K+1)-1) ;
hess.submat(dd*(K+1), (nD+1)*(K+1), (dd+1)*(K+1)-1, (nD+1)*(K+1)+numProdsTotal-1) = trans(dL2ddeltadtheta2.cols((dd-1)*(K+1), dd*(K+1)-1)) ;
} 
 
//HESSIAN FOR DELTA 
hess.submat(((nD+1)*(K+1)), ((nD+1)*(K+1)), ((nD+1)*(K+1)+numProdsTotal-1), ((nD+1)*(K+1)+numProdsTotal-1)) = dL2ddelta2 ;

return hess ;
} 
  

// [[Rcpp::export]]
List blk_cpp( mat simShare, mat multip, int nn, mat dL2ddelta2DIAG_old, mat dL2ddeltadtheta_old, mat dL2dtheta2_old, int index_start, int index_end, mat x_char_index, mat x_char, mat ooo1, mat ooo, mat v, int K, int nD, int prodsMarket_t ) {
  
mat dL2ddelta2DIAG_new = dL2ddelta2DIAG_old; 
mat dL2ddeltadtheta2_new = dL2ddeltadtheta_old;
mat dL2dtheta2_new = dL2dtheta2_old;

 
for ( int rr = 1; rr < nn+1; ++rr ) { 
  
mat simS_mat = simShare.submat(index_start-1, rr-1, index_end-1, rr-1);  
int rows = simS_mat.n_rows;  
mat simS_diag = eye(rows,rows)%repmat(simS_mat,1,rows);
double sumprod_mpsimS = as_scalar(trans(multip)*simS_mat);
mat v_rr = v.cols(rr-1,rr-1);
mat xsimSx = x_char_index - (ooo1*(trans(simS_mat)*x_char_index));

mat xsimSxv_mat = zeros(rows,15) ;
mat dSdtheta2rr_mat = zeros(rows,15);
for ( int dd = 1; dd < (nD+1); ++dd ) {       
xsimSxv_mat.cols((K+1)*(dd-1), (K+1)*dd-1) = xsimSx%(ooo1*trans(v.submat((K+1)*(dd-1), rr-1, dd*(K+1)-1, rr-1))) ;  
dSdtheta2rr_mat.cols((K+1)*(dd-1), (K+1)*dd-1) = (simS_mat*ooo)%xsimSxv_mat.cols((K+1)*(dd-1), (K+1)*dd-1) ;
}


mat blk12 = sumprod_mpsimS*(-simS_diag + 2*simS_mat*trans(simS_mat)) -(multip%simS_mat)*trans(simS_mat) - trans((multip%simS_mat)*trans(simS_mat));      
int rows_multip = multip.n_rows;  
mat blk3 = eye(rows_multip,rows_multip)%repmat(multip%simS_mat,1,rows_multip);  
mat blk = blk12 + blk3; 


dL2ddelta2DIAG_new.submat(index_start-1, 0, index_end-1, prodsMarket_t-1) += blk/nn;
mat dL2ddeltadtheta2rr_dd ;
mat dL2dtheta2rr_dd ;
for ( int dd = 1; dd < (nD+1); ++dd ) { 
dL2ddeltadtheta2rr_dd = -simS_mat*trans(multip)*dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1) - sumprod_mpsimS*dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1) + (multip*ooo)%dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1) ;   
dL2dtheta2rr_dd = trans(((multip*ooo)%dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1)))*xsimSxv_mat.cols((K+1)*(dd-1), dd*(K+1)-1) + sumprod_mpsimS*(-trans(dSdtheta2rr_mat.cols((K+1)*(dd-1), dd*(K+1)-1))*x_char.rows(index_start-1, index_end-1)%(ones(K+1,1)*trans(v.submat((K+1)*(dd-1), rr-1, dd*(K+1)-1, rr-1)))) ;
dL2ddeltadtheta2_new.submat(index_start-1,((dd-1)*(K+1)), index_end-1, (dd*(K+1)-1)) += dL2ddeltadtheta2rr_dd/nn;
dL2dtheta2_new.submat(((dd-1)*(K+1)), ((dd-1)*(K+1)), (dd*(K+1))-1, (dd*(K+1))-1) += dL2dtheta2rr_dd/nn;
mat dL2dtheta2rr_dd_ddd;
for ( int ddd = 1; ddd < (nD+1); ++ddd ) {  
if (ddd!=dd) {   
dL2dtheta2rr_dd_ddd = trans((multip*ooo)%dSdtheta2rr_mat.cols((K+1)*(ddd-1), ddd*(K+1)-1))*xsimSxv_mat.cols((K+1)*(dd-1), dd*(K+1)-1) + sumprod_mpsimS*(-trans(dSdtheta2rr_mat.cols((K+1)*(ddd-1), ddd*(K+1)-1))*x_char.rows(index_start-1, index_end-1)%(ones(K+1,1)*trans(v.submat((K+1)*(dd-1), rr-1, dd*(K+1)-1, rr-1)))) ;   
dL2dtheta2_new.submat(((dd-1)*(K+1)), ((ddd-1)*(K+1)), (dd*(K+1))-1, (ddd*(K+1))-1) += dL2dtheta2rr_dd_ddd/nn ;
}
}
}


}

return List::create( Named("dL2ddelta2DIAG_new") = dL2ddelta2DIAG_new, Named("dL2ddeltadtheta2_new") = dL2ddeltadtheta2_new, Named("dL2dtheta2_new") = dL2dtheta2_new ) ; 
} 



// [[Rcpp::export]]
List blk_cpp_new( mat simShare, mat multip, int nn, mat dL2ddelta2DIAG_old, mat dL2ddeltadtheta_old, mat dL2dtheta2_old, int index_start, int index_end, mat x_char_index, mat x_char, mat ooo1, mat ooo, mat v, int K, int nD, int prodsMarket_t, int t ) {
  
mat dL2ddelta2DIAG_new = dL2ddelta2DIAG_old; 
mat dL2ddeltadtheta2_new = dL2ddeltadtheta_old;
mat dL2dtheta2_new = dL2dtheta2_old;

 
for ( int rr = 1; rr < nn+1; ++rr ) { 
  
mat simS_mat = simShare.submat(index_start-1, rr-1, index_end-1, rr-1);  
int rows = simS_mat.n_rows;  
mat simS_diag = eye(rows,rows)%repmat(simS_mat,1,rows);
double sumprod_mpsimS = as_scalar(trans(multip)*simS_mat);
mat v_rr = v.cols(rr-1,rr-1);
mat xsimSx = x_char_index - (ooo1*(trans(simS_mat)*x_char_index));


// v    
mat xsimSxv_v = xsimSx%(ooo1*trans(v.submat(0, rr-1, (K+1)-1, rr-1))) ;  
mat dSdtheta2rr_v = (simS_mat*ooo)%xsimSxv_v ;
// Y
mat xsimSxv_Y = xsimSx%(ooo1*trans(v.submat((K+1), rr-1, 2*(K+1)-1, rr-1))) ;  
mat dSdtheta2rr_Y = (simS_mat*ooo)%xsimSxv_Y ;


mat blk12 = sumprod_mpsimS*(-simS_diag + 2*simS_mat*trans(simS_mat)) -(multip%simS_mat)*trans(simS_mat) - trans((multip%simS_mat)*trans(simS_mat));      
int rows_multip = multip.n_rows;  
mat blk3 = eye(rows_multip,rows_multip)%repmat(multip%simS_mat,1,rows_multip);  
mat blk = blk12 + blk3; 


dL2ddelta2DIAG_new.submat(index_start-1, 0, index_end-1, prodsMarket_t-1) += blk/nn;
mat dL2ddeltadtheta2rr ;
mat dL2dtheta2rr ;
mat dL2dtheta2rr_v ;
mat dL2dtheta2rr_Y ;


dL2ddeltadtheta2rr = -simS_mat*trans(multip)*dSdtheta2rr_v - sumprod_mpsimS*dSdtheta2rr_v + (multip*ooo)%dSdtheta2rr_v ; 
dL2ddeltadtheta2_new.submat(index_start-1, 0, index_end-1, (K+1)-1) += dL2ddeltadtheta2rr/nn;
dL2ddeltadtheta2rr = -simS_mat*trans(multip)*dSdtheta2rr_Y - sumprod_mpsimS*dSdtheta2rr_Y + (multip*ooo)%dSdtheta2rr_Y ; 
dL2ddeltadtheta2_new.submat(index_start-1, (K+1), index_end-1, 2*(K+1)-1) += dL2ddeltadtheta2rr/nn;

dL2dtheta2rr_v = trans(((multip*ooo)%dSdtheta2rr_v))*xsimSxv_v + sumprod_mpsimS*(-trans(dSdtheta2rr_v)*x_char.rows(index_start-1, index_end-1)%(ones(K+1,1)*trans(v.submat(0, rr-1, (K+1)-1, rr-1)))) ;
dL2dtheta2rr_Y = trans(((multip*ooo)%dSdtheta2rr_Y))*xsimSxv_Y + sumprod_mpsimS*(-trans(dSdtheta2rr_Y)*x_char.rows(index_start-1, index_end-1)%(ones(K+1,1)*trans(v.submat((K+1), rr-1, 2*(K+1)-1, rr-1)))) ; 
dL2dtheta2_new.submat(0, 0, (K+1)-1, (K+1)-1) += dL2dtheta2rr_v/nn;
dL2dtheta2_new.submat(K+1, K+1, 2*(K+1)-1, 2*(K+1)-1) += dL2dtheta2rr_Y/nn;

dL2dtheta2rr = trans((multip*ooo)%dSdtheta2rr_Y)*xsimSxv_v + sumprod_mpsimS*(-trans(dSdtheta2rr_Y)*x_char.rows(index_start-1, index_end-1)%(ones(K+1,1)*trans(v.submat(0, rr-1, (K+1)-1, rr-1)))) ;   
dL2dtheta2_new.submat(0, (K+1), (K+1)-1, 2*(K+1)-1) += dL2dtheta2rr/nn ;

dL2dtheta2rr = trans((multip*ooo)%dSdtheta2rr_v)*xsimSxv_Y + sumprod_mpsimS*(-trans(dSdtheta2rr_v)*x_char.rows(index_start-1, index_end-1)%(ones(K+1,1)*trans(v.submat((K+1), rr-1, 2*(K+1)-1, rr-1)))) ; 
dL2dtheta2_new.submat(K+1, 0, 2*(K+1)-1, (K+1)-1) += dL2dtheta2rr/nn ;
 
}

return List::create( Named("dL2ddelta2DIAG_new") = dL2ddelta2DIAG_new, Named("dL2ddeltadtheta2_new") = dL2ddeltadtheta2_new, Named("dL2dtheta2_new") = dL2dtheta2_new ) ; 
} 



// JACOBIAN
// [[Rcpp::export]]
mat dc_1_cpp(mat simShare, int index_start, int index_end, mat dSddeltaDIAG_old, int nn ) {
 
mat dSddeltaDIAG_temp = dSddeltaDIAG_old; 
 
for ( int rr = 1; rr < nn+1; ++rr ) {    
    
mat simShare_index_rr = simShare.submat(index_start-1, rr-1, index_end-1, rr-1);
int rows = simShare_index_rr.n_rows; 
mat diag_simShare_index_rr = eye(rows, rows)%repmat(simShare_index_rr,1,rows);


dSddeltaDIAG_temp += (diag_simShare_index_rr - simShare_index_rr*trans(simShare_index_rr))/nn; 
}

return dSddeltaDIAG_temp ;
} 


// [[Rcpp::export]]
mat dc_2_cpp(mat dSdtheta2_old, mat simShare, mat ooo, mat ooo1, mat v, mat x_char, int nn, int index_start, int index_end, int nD, int K) {
        
mat dSdtheta2_new = dSdtheta2_old ;

for ( int rr = 1; rr < nn+1; ++rr ) {
for ( int dd = 1; dd < (nD+1); ++dd ) { 
dSdtheta2_new.submat(index_start-1, (dd-1)*(K+1), index_end-1, dd*(K+1)-1) += (simShare.submat(index_start-1, rr-1, index_end-1, rr-1)*ooo)%(ooo1*trans(v.submat((K+1)*(dd-1), rr-1, dd*(K+1)-1, rr-1)))%( x_char.rows(index_start-1,index_end-1) - (ooo1*(trans(simShare.submat(index_start-1, rr-1, index_end-1, rr-1))*x_char.rows(index_start-1,index_end-1))))/nn ;
}
}

return dSdtheta2_new ;
} 


// [[Rcpp::export]]
mat dc_2_cpp_new(mat dSdtheta2_old, mat simShare, mat ooo, mat ooo1, mat v, mat x_char, int nn, int index_start, int index_end, int nD, int K, int K_A, int t) {
        
mat dSdtheta2_new = dSdtheta2_old ;

for ( int rr = 1; rr < nn+1; ++rr ) {
for ( int dd = 1; dd < nD+1; ++dd ) {   
dSdtheta2_new.submat(index_start-1, ((dd-1)*(K_A+1)), index_end-1, dd*(K_A+1)-1) += (simShare.submat(index_start-1, rr-1, index_end-1, rr-1)*ooo)%(ooo1*trans(v.submat(((dd-1)*(K+1)), rr-1, ((dd-1)*(K+1)+K_A), rr-1)))%( x_char.submat(index_start-1, 0, index_end-1, K_A) - (ooo1*(trans(simShare.submat(index_start-1, rr-1, index_end-1, rr-1))*x_char.submat(index_start-1, 0, index_end-1, K_A))))/nn ;
}
}

return dSdtheta2_new ;
} 


// [[Rcpp::export]]
mat dSdbetaf_cpp(mat dSdbetaf_old, mat sigmaf, mat simShare, mat ooo, mat ooo1, mat v, mat x_char, int nn, int index_start, int index_end, int nD, int K, int K_A, int t, int countries) {
       
mat dSdbetaf_new = dSdbetaf_old ;  
 
for ( int rr = 1; rr < nn+1; ++rr ) {
mat simShare_index_rr = simShare.submat(index_start-1, rr-1, index_end-1, rr-1);
int rows = simShare_index_rr.n_rows; 
mat diag_simShare_index_rr = eye(rows, rows)%repmat(simShare_index_rr,1,rows); 
dSdbetaf_new += ((diag_simShare_index_rr - simShare_index_rr*trans(simShare_index_rr)) * x_char.submat(index_start-1, K_A+1, index_end-1, K))/nn;  
for ( int dd = 1; dd < nD+1; ++dd ) {   
mat v_rr = v.submat(((dd-1)*(K+1)+ (K_A+1)), rr-1, dd*(K+1)-1,rr-1);
mat diag_v_rr = eye(v_rr.n_rows, v_rr.n_rows)%repmat(v_rr,1,v_rr.n_rows);   

mat sigmaf_dd = repmat(sigmaf.rows(dd-1,dd-1), v_rr.n_rows, v_rr.n_rows);
dSdbetaf_new += ((diag_simShare_index_rr - simShare_index_rr*trans(simShare_index_rr)) * x_char.submat(index_start-1, K_A+1, index_end-1, K) * (diag_v_rr % sigmaf_dd))/nn;   
}   
}
        
return dSdbetaf_new ;
} 







// [[Rcpp::export]]
mat dSdsigmaf_cpp(mat dSdsigmaf_old, mat betaf, mat sigmaf, mat simShare, mat ooo, mat ooo1, mat ooo2, mat v, mat x_char, int nn, int index_start, int index_end, int nD, int K, int K_A, int t, int countries) {
       
mat dSdsigmaf_new = dSdsigmaf_old ;  
 
for ( int rr = 1; rr < nn+1; ++rr ) {

mat vf_mat = v.submat(K_A+1, rr-1, K, rr-1);
for ( int dd = 2; dd < nD+1; ++dd ) {
vf_mat = join_rows(vf_mat, v.submat((dd-1)*(K+1) + (K_A+1), rr-1, dd*(K+1)-1, rr-1));
}
mat simShare_index_rr = simShare.submat(index_start-1, rr-1, index_end-1, rr-1);
dSdsigmaf_new += ((simShare_index_rr*ooo2) % ((x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)*vf_mat) % (x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)*betaf*ooo2) - ((ooo1*trans(simShare_index_rr)) * (x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)*betaf*ooo2) % (x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)*vf_mat))))/nn;     
  
}
        
return dSdsigmaf_new ;

} 




// [[Rcpp::export]]
List blk_cpp_FE( mat simShare, mat multip, int nn, mat dL2ddelta2DIAG_old, mat dL2ddeltadtheta_old, mat dL2dtheta2_old, mat dL2ddeltadbetaf_old, mat dL2dbetaf2_old, mat dL2ddeltadsigmaf_old, mat dL2dbetafdtheta2_old, mat dL2dbetafdsigmaf_old, mat dL2dsigmaf2_old, mat dL2dtheta2dsigmaf_old, int index_start, int index_end, mat x_char_index, mat x_char, mat betaf, mat sigmaf, mat v_market_v, mat v_market_Y, mat b_market, mat ooo1, mat ooo3, mat v, int K, int K_A, int nD, int prodsMarket_t, int t ) {

mat ooo = ones(1,K_A+1);
mat ooo2 = ones(1,nD);


mat dL2ddelta2DIAG_new = dL2ddelta2DIAG_old; 
mat dL2ddeltadtheta2_new = dL2ddeltadtheta_old;
mat dL2dtheta2_new = dL2dtheta2_old;

mat dL2ddeltadbetaf_new = dL2ddeltadbetaf_old;
mat dL2dbetaf2_new = dL2dbetaf2_old;
mat dL2ddeltadsigmaf_new = dL2ddeltadsigmaf_old;
mat dL2dbetafdtheta2_new = dL2dbetafdtheta2_old;
mat dL2dbetafdsigmaf_new = dL2dbetafdsigmaf_old;
mat dL2dsigmaf2_new = dL2dsigmaf2_old;
mat dL2dtheta2dsigmaf_new = dL2dtheta2dsigmaf_old;
 
for ( int rr = 1; rr < nn+1; ++rr ) { 
  
mat simS_mat = simShare.submat(index_start-1, rr-1, index_end-1, rr-1);
mat vf_mat = v.submat(K_A+1, rr-1, K, rr-1);
for ( int dd = 2; dd < nD+1; ++dd ) {
vf_mat = join_rows(vf_mat, v.submat((dd-1)*(K+1) + (K_A+1), rr-1, dd*(K+1)-1, rr-1));
}
int rows = simS_mat.n_rows;  
mat simS_diag = eye(rows,rows)%repmat(simS_mat,1,rows);
double sumprod_mpsimS = as_scalar(trans(multip)*simS_mat);
mat v_rr = v.cols(rr-1,rr-1);
mat xsimSx = x_char.submat(index_start-1, 0, index_end-1, K_A) - (ooo1*(trans(simS_mat)*x_char.submat(index_start-1, 0, index_end-1, K_A)));


// v    
mat xsimSxv_v = xsimSx%(ooo1*trans(v.submat(0, rr-1, (K_A+1)-1, rr-1))) ;  
mat bsimSbv_v = b_market%v_market_v.cols(rr-1,rr-1) - ooo1*trans(simS_mat)*(b_market%v_market_v.cols(rr-1, rr-1));
mat dSdtheta2rr_v = (simS_mat*ooo)%xsimSxv_v ;
// Y
mat xsimSxv_Y = xsimSx%(ooo1*trans(v.submat((K+1), rr-1, (K+1)+K_A, rr-1))) ;  
mat bsimSbv_Y = b_market%v_market_Y.cols(rr-1,rr-1) - ooo1*trans(simS_mat)*(b_market%v_market_Y.cols(rr-1, rr-1));
mat dSdtheta2rr_Y = (simS_mat*ooo)%xsimSxv_Y ;


mat blk12 = sumprod_mpsimS*(-simS_diag + 2*simS_mat*trans(simS_mat)) -(multip%simS_mat)*trans(simS_mat) - trans((multip%simS_mat)*trans(simS_mat));      
int rows_multip = multip.n_rows;  
mat blk3 = eye(rows_multip,rows_multip)%repmat(multip%simS_mat,1,rows_multip);  
mat blk = blk12 + blk3; 


dL2ddelta2DIAG_new.submat(index_start-1, 0, index_end-1, prodsMarket_t-1) += blk/nn;
dL2ddeltadbetaf_new.rows(index_start-1, index_end-1) += (trans((vf_mat*sigmaf + ooo3)*trans(ooo1))%(blk*x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)))/nn;
dL2dbetaf2_new += (((vf_mat*sigmaf + ooo3)*trans(ooo3)) % ((vf_mat*sigmaf + ooo3)*trans(ooo3)) % (trans(x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1))*blk*x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)))/nn;

mat dSdsigmafrr = (simS_mat*ooo2) % ((x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)*vf_mat) % (x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)*betaf*ooo2) - ((ooo1*trans(simS_mat)) * ((x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)*betaf*ooo2) % (x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1)*vf_mat))));
        
mat dL2ddeltadtheta2rr;
mat dL2ddeltadsigmafrr;
mat dL2dbetafdtheta2rr;
mat dL2dbetafdsigmafrr;
mat dL2dtheta2rr;
mat dL2dtheta2rr_v;
mat dL2dtheta2rr_Y;
mat dL2dsigmaf2rr;
mat dL2dtheta2dsigmafrr;

dL2ddeltadtheta2rr = -simS_mat*trans(multip)*dSdtheta2rr_v - sumprod_mpsimS*dSdtheta2rr_v + (multip*ooo)%dSdtheta2rr_v ; 
dL2dbetafdtheta2rr = ((vf_mat*sigmaf + ooo3)*ooo) % (trans(x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1))*dL2ddeltadtheta2rr);
dL2ddeltadtheta2_new.submat(index_start-1, 0, index_end-1, (K_A+1)-1) += dL2ddeltadtheta2rr/nn;
dL2dbetafdtheta2_new.cols(0, K_A) += dL2dbetafdtheta2rr/nn;
dL2ddeltadtheta2rr = -simS_mat*trans(multip)*dSdtheta2rr_Y - sumprod_mpsimS*dSdtheta2rr_Y + (multip*ooo)%dSdtheta2rr_Y ; 
dL2dbetafdtheta2rr = ((vf_mat*sigmaf + ooo3)*ooo) % (trans(x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1))*dL2ddeltadtheta2rr);
dL2ddeltadtheta2_new.submat(index_start-1, (K_A+1), index_end-1, 2*(K_A+1)-1) += dL2ddeltadtheta2rr/nn;
dL2dbetafdtheta2_new.cols(K_A+1, 2*(K_A+1)-1) += dL2dbetafdtheta2rr/nn;



for ( int dd = 1; dd < nD+1; ++dd ) {
dL2ddeltadsigmafrr = -simS_mat*trans(multip)*dSdsigmafrr.cols(dd-1,dd-1) - sumprod_mpsimS*dSdsigmafrr.cols(dd-1,dd-1) + multip%dSdsigmafrr.cols(dd-1,dd-1);
dL2dbetafdsigmafrr = (trans(x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1))*dL2ddeltadsigmafrr) % (vf_mat*sigmaf + ooo3) + (trans(x_char.submat(index_start-1, K_A+1, index_end-1, x_char.n_cols-1))*(simS_mat%(multip - sumprod_mpsimS*ooo1))) % vf_mat.cols(dd-1,dd-1);
dL2ddeltadsigmaf_new.submat(index_start-1, dd-1, index_end-1, dd-1) += dL2ddeltadsigmafrr/nn;
dL2dbetafdsigmaf_new.cols(dd-1, dd-1) += dL2dbetafdsigmafrr/nn; 
}


dL2dtheta2rr_v = trans(((multip*ooo)%dSdtheta2rr_v))*xsimSxv_v + sumprod_mpsimS*(-trans(dSdtheta2rr_v)*x_char.submat(index_start-1, 0, index_end-1, K_A)%(ones(K_A+1,1)*trans(v.submat(0, rr-1, (K_A+1)-1, rr-1)))) ;
dL2dtheta2rr_Y = trans(((multip*ooo)%dSdtheta2rr_Y))*xsimSxv_Y + sumprod_mpsimS*(-trans(dSdtheta2rr_Y)*x_char.submat(index_start-1, 0, index_end-1, K_A)%(ones(K_A+1,1)*trans(v.submat(K+1, rr-1, (K+1) + (K_A+1)-1, rr-1)))) ;
dL2dtheta2_new.submat(0, 0, K_A, K_A) += dL2dtheta2rr_v/nn;
dL2dtheta2_new.submat(K_A+1, K_A+1, 2*(K_A+1)-1, 2*(K_A+1)-1) += dL2dtheta2rr_Y/nn;


dL2dsigmaf2rr = trans(multip%dSdsigmafrr.cols(0,0))*bsimSbv_v + sumprod_mpsimS*(-trans(dSdsigmafrr.cols(0,0))*(b_market%v_market_v.cols(rr-1,rr-1)));
dL2dsigmaf2_new.submat(0,0,0,0) += dL2dsigmaf2rr/nn;
dL2dsigmaf2rr = trans(multip%dSdsigmafrr.cols(1,1))*bsimSbv_Y + sumprod_mpsimS*(-trans(dSdsigmafrr.cols(1,1))*(b_market%v_market_Y.cols(rr-1,rr-1)));
dL2dsigmaf2_new.submat(1,1,1,1) += dL2dsigmaf2rr/nn;

dL2dtheta2rr = trans((multip*ooo)%dSdtheta2rr_Y)*xsimSxv_v + sumprod_mpsimS*(-trans(dSdtheta2rr_Y)*x_char.submat(index_start-1, 0, index_end-1, K_A)%(ones(K_A+1,1)*trans(v.submat(0, rr-1, (K_A+1)-1, rr-1)))) ;   
dL2dtheta2_new.submat(0, K_A+1, K_A, 2*(K_A+1)-1) += dL2dtheta2rr/nn ;

dL2dtheta2rr = trans((multip*ooo)%dSdtheta2rr_v)*xsimSxv_Y + sumprod_mpsimS*(-trans(dSdtheta2rr_v)*x_char.submat(index_start-1, 0, index_end-1, K_A)%(ones(K_A+1,1)*trans(v.submat(K+1, rr-1, (K+1) + (K_A+1)-1, rr-1)))) ; 
dL2dtheta2_new.submat(K_A+1, 0, 2*(K_A+1)-1, K_A) += dL2dtheta2rr/nn ;
 
 
dL2dsigmaf2rr = trans(multip%dSdsigmafrr.cols(1,1))*bsimSbv_v + sumprod_mpsimS*(-trans(dSdsigmafrr.cols(1,1))*(b_market%v_market_v.cols(rr-1,rr-1)));
dL2dsigmaf2_new.submat(0,1,0,1) += dL2dsigmaf2rr/nn;
dL2dsigmaf2rr = trans(multip%dSdsigmafrr.cols(0,0))*bsimSbv_Y + sumprod_mpsimS*(-trans(dSdsigmafrr.cols(0,0))*(b_market%v_market_Y.cols(rr-1,rr-1)));
dL2dsigmaf2_new.submat(1,0,1,0) += dL2dsigmaf2rr/nn;


//dd=0 ddd=0
dL2dtheta2dsigmafrr = trans(multip%dSdsigmafrr.cols(0,0))*xsimSxv_v + sumprod_mpsimS*(-trans(dSdsigmafrr.cols(0,0))*x_char.submat(index_start-1, 0, index_end-1, K_A)%(trans(v.submat(0, rr-1, K_A, rr-1)))) ;    
dL2dtheta2dsigmaf_new.submat(0, 0, 0, K_A) += dL2dtheta2dsigmafrr/nn;
//dd=1 ddd=1
dL2dtheta2dsigmafrr = trans(multip%dSdsigmafrr.cols(1,1))*xsimSxv_Y + sumprod_mpsimS*(-trans(dSdsigmafrr.cols(1,1))*x_char.submat(index_start-1, 0, index_end-1, K_A)%(trans(v.submat(K+1, rr-1, (K+1) + K_A, rr-1)))) ;    
dL2dtheta2dsigmaf_new.submat(1, K_A+1, 1, 2*(K_A+1)-1) += dL2dtheta2dsigmafrr/nn;
//dd=1 ddd=0
dL2dtheta2dsigmafrr = trans(multip%dSdsigmafrr.cols(0,0))*xsimSxv_Y + sumprod_mpsimS*(-trans(dSdsigmafrr.cols(0,0))*x_char.submat(index_start-1, 0, index_end-1, K_A)%(trans(v.submat(K+1, rr-1, (K+1) + K_A, rr-1)))) ;    
dL2dtheta2dsigmaf_new.submat(0, K_A+1, 0, 2*(K_A+1)-1) += dL2dtheta2dsigmafrr/nn;
//dd=0 ddd=1
dL2dtheta2dsigmafrr = trans(multip%dSdsigmafrr.cols(1,1))*xsimSxv_v + sumprod_mpsimS*(-trans(dSdsigmafrr.cols(1,1))*x_char.submat(index_start-1, 0, index_end-1, K_A)%(trans(v.submat(0, rr-1, K_A, rr-1)))) ;    
dL2dtheta2dsigmaf_new.submat(1, 0, 1, K_A) += dL2dtheta2dsigmafrr/nn;

}

return List::create( Named("dL2ddelta2DIAG_new") = dL2ddelta2DIAG_new, Named("dL2ddeltadtheta2_new") = dL2ddeltadtheta2_new, Named("dL2dtheta2_new") = dL2dtheta2_new, Named("dL2ddeltadbetaf_new") = dL2ddeltadbetaf_new, Named("dL2dbetaf2_new") = dL2dbetaf2_new, Named("dL2ddeltadsigmaf_new") = dL2ddeltadsigmaf_new, Named("dL2dbetafdtheta2_new") = dL2dbetafdtheta2_new, Named("dL2dbetafdsigmaf_new") = dL2dbetafdsigmaf_new, Named("dL2dsigmaf2_new") = dL2dsigmaf2_new, Named("dL2dtheta2dsigmaf_new") = dL2dtheta2dsigmaf_new) ; 
} 



// [[Rcpp::export]]
mat hess_cpp_FE(mat hess_old, mat dL2dtheta2, mat dL2ddeltadtheta2, mat dL2ddelta2, mat dL2dbetaf2, mat dL2dbetafdtheta2, mat dL2dbetafdsigmaf, mat dL2ddeltadbetaf, mat dL2dsigmaf2, mat dL2dtheta2dsigmaf, mat dL2ddeltadsigmaf, int nD, int K, int K_A, int countries, int numProdsTotal) { 
  
mat hess = hess_old ;  
 
// HESSIAN FOR THETA2, THETA3, THETA4 
for ( int dd = 1; dd < (nD+1); ++dd ) {
hess.submat(dd*(K_A+1), dd*(K_A+1), (dd+1)*(K_A+1)-1, (dd+1)*(K_A+1)-1) = dL2dtheta2.submat((dd-1)*(K_A+1), (dd-1)*(K_A+1), dd*(K_A+1)-1, dd*(K_A+1)-1) ;    
} 
 
for ( int dd = 1; dd < (nD+1); ++dd ) {  
for ( int ddd = 1; ddd < (nD+1); ++ddd ) {  
if (ddd!=dd) {     
hess.submat(dd*(K_A+1), ddd*(K_A+1), (dd+1)*(K_A+1)-1, (ddd+1)*(K_A+1)-1) = trans(dL2dtheta2.submat((dd-1)*(K_A+1), (ddd-1)*(K_A+1), dd*(K_A+1)-1, ddd*(K_A+1)-1)) ;    
}
}
} 

//Hessian for betaf
hess.submat((nD+1)*(K_A+1), (nD+1)*(K_A+1), (nD+1)*(K_A+1)+(countries-1)-1, (nD+1)*(K_A+1)+(countries-1)-1) = dL2dbetaf2;
hess.submat((nD+1)*(K_A+1), K_A+1, (nD+1)*(K_A+1)+(countries-1)-1, (nD+1)*(K_A+1)-1) = dL2dbetafdtheta2;    
hess.submat(K_A+1, (nD+1)*(K_A+1), (nD+1)*(K_A+1)-1, (nD+1)*(K_A+1)+(countries-1)-1) = trans(dL2dbetafdtheta2);   
hess.submat((nD+1)*(K_A+1), (nD+1)*(K_A+1)+(countries-1), (nD+1)*(K_A+1)+(countries-1)-1, (nD+1)*(K_A+1)+(countries-1)+nD-1) = dL2dbetafdsigmaf; 
hess.submat((nD+1)*(K_A+1)+(countries-1), (nD+1)*(K_A+1), (nD+1)*(K_A+1)+(countries-1)+nD-1, (nD+1)*(K_A+1)+(countries-1)-1) = trans(dL2dbetafdsigmaf); 
hess.submat((nD+1)*(K_A+1), (nD+1)*(K_A+1)+(countries-1)+nD, (nD+1)*(K_A+1)+(countries-1)-1, (nD+1)*(K_A+1)+(countries-1)+nD + numProdsTotal - 1) = trans(dL2ddeltadbetaf);  
hess.submat((nD+1)*(K_A+1)+(countries-1)+nD, (nD+1)*(K_A+1), (nD+1)*(K_A+1)+(countries-1)+nD + numProdsTotal - 1, (nD+1)*(K_A+1)+(countries-1)-1) = dL2ddeltadbetaf; 
 
//Hessian for sigmaf
hess.submat((nD+1)*(K_A+1)+(countries-1), (nD+1)*(K_A+1)+(countries-1), (nD+1)*(K_A+1)+(countries-1)+nD-1, (nD+1)*(K_A+1)+(countries-1)+nD-1) = dL2dsigmaf2;    
hess.submat((nD+1)*(K_A+1)+(countries-1), K_A+1, (nD+1)*(K_A+1)+(countries-1)+nD-1, (nD+1)*(K_A+1)-1) = dL2dtheta2dsigmaf;  
hess.submat(K_A+1, (nD+1)*(K_A+1)+(countries-1), (nD+1)*(K_A+1)-1, (nD+1)*(K_A+1)+(countries-1)+nD-1) = trans(dL2dtheta2dsigmaf); 
hess.submat((nD+1)*(K_A+1)+(countries-1), (nD+1)*(K_A+1)+(countries-1)+nD, (nD+1)*(K_A+1)+(countries-1)+nD-1, (nD+1)*(K_A+1)+(countries-1)+nD + numProdsTotal-1) = trans(dL2ddeltadsigmaf); 
hess.submat((nD+1)*(K_A+1)+(countries-1)+nD, (nD+1)*(K_A+1)+(countries-1), (nD+1)*(K_A+1)+(countries-1)+nD + numProdsTotal-1, (nD+1)*(K_A+1)+(countries-1)+nD-1) = dL2ddeltadsigmaf;


// HESSIAN FOR DELTA AND THETA2dd
for ( int dd = 1; dd < (nD+1); ++dd ) {
hess.submat((nD+1)*(K_A+1) + (countries-1) + nD, dd*(K_A+1), (nD+1)*(K_A+1) + (countries-1) + nD + numProdsTotal-1, (dd+1)*(K_A+1) -1) = dL2ddeltadtheta2.cols((dd-1)*(K_A+1), dd*(K_A+1)-1);
hess.submat(dd*(K_A+1), (nD+1)*(K_A+1) + (countries-1) + nD, (dd+1)*(K_A+1)-1, (nD+1)*(K_A+1) + (countries-1) + nD + numProdsTotal -1) = trans(dL2ddeltadtheta2.cols((dd-1)*(K_A+1), dd*(K_A+1)-1));
} 
 
//HESSIAN FOR DELTA 
hess.submat((nD+1)*(K_A+1) + (countries-1) + nD, (nD+1)*(K_A+1) + (countries-1) + nD, (nD+1)*(K_A+1) + (countries-1) + nD + numProdsTotal - 1, (nD+1)*(K_A+1) + (countries-1) + nD + numProdsTotal -1) = dL2ddelta2; 

return hess ;
} 


// [[Rcpp::export]]
mat zeros_test(int N, int M) {
       
mat gh = zeros(M,N) ;  
  
return gh ;

} 
